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Abstract - Polynomial Chaos Expansions represent a 
powerful tool to simulate stochastic models of dynam- 
ical systems. Yet, deriving the expansion's coefficients 
for complex systems might require a significant and non- 
trivial manipulation of the model, or the computation 
of large numbers of simulation runs, rendering the ap- 
proach too time consuming and impracticable for ap- 
plications with more than a handful of random vari- 
ables. We introduce a novel computationally tractable 
technique for computing the coefficients of polynomial 
chaos expansions. The approach exploits a regulariza- 
tion technique with a particular choice of weighting ma- 
trices, which allow to take into account the specific fea- 
tures of Polynomial Chaos expansions. The method, com- 
pletely based on convex optimization, can be applied to 
problems with a large number of random variables and 
uses a modest number of Monte Carlo simulations, while 
avoiding model manipulations. Additional information 
on the stochastic process, when available, can be also 
incorporated in the approach by means of convex con- 
straints. We show the effectiveness of the proposed tech- 
nique in three appUcations in diverse fields, including 
the analysis of a nonUnear electric circuit, a chaotic mo- 
del of organizational behavior, finally a chemical oscilla- 
tor. 
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1 Introduction 

In most science and engineering applications, there is the 
need to simulate mathematical models of the process under 
study, in the form of ordinary or partial differential equa- 
tions, with the aim to figure out the time (or space) course 
of a variable of interest v for analysis, decision-making and 
control. In many cases, in order to take into account the 
presence of uncertainty, unknown external inputs and in 
general any effect that produces a mismatch between the 
model equations and reality, the model to be simulated is 
not fully deterministic: uncertainty and disturbance are of- 
ten modeled as quantities with stochastic nature, named here 
"input random variables" and indicated with 6, hence the 
name stochastic models. In these cases, v{0) is a random 
variable, too, and one is interested in computing its statis- 
tics. The issue of simulating complex, nonlinear stochas- 
tic models with sufficiently high accuracy and low compu- 
tational effort is still a challenge in important and diverse 
fields, like analysis of large power grids, weather forecasts 
at different scales and simulation of biological systems, to 
name just a few. The typical approach followed to simulate 
a stochastic model is the well-known Monte Carlo (MC) 
technique, which relies on the sampling of a finite number 
M of values of 0, according to its distribution. With suf- 
ficiently large M, the MC approach gives good statistical 
estimates (e.g. first and second order moments) of the vari- 
ables of interest, and also of its probability density function 
(pdf). However, the application of MC simulations with 
the system model may be too computationally demanding, 
particularly in those cases when the model is complex and 
the inherent variables have large dimensions. Polynomial 
Chaos Expansions (PCEs) (see e.g. |[l]|2l[3l|4l|5]|6l|2l|8ll3) 
provide a useful tool to significantly reduce the computa- 
tional effort required to simulate a stochastic system, by 
conceptually replacing the mapping between 6 and v, im- 
plicitly defined by the integration of the model's differen- 
tial equations, with an explicit function v{6), which takes 
the form of a truncated series of polynomials. The poly- 
nomials in the PCE are orthogonal, so that the statistical 
moments of v{d) can be computed directly from the ex- 
pansion's coefficients. Moreover, the computational effort 
required to evaluate a PCE is often orders of magnitude 
lower than the one required to simulate the system model: 
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therefore, it is possible to estimate the pdf of v by using 
a Monte Carlo approach with the PCE, instead of the sys- 
tem model, with significant time savings. As an example, in 
the first case study considered in this paper, 100,000 Monte 
Carlo simulations with the system model require 3,800 s, 
while the same number of PCE evaluations are obtained 
in 24 s, with the same hardware and software. PCEs have 
been used with good results in a number of different ar- 
eas, including experimental modeling, materials, mechan- 
ics, power systems, systems biology, and control, see e.g. 
||T0l[IT][l3|T3[ll|Bl|I3[I3 18 19 2Q]|n]l22l. Yet, 
the computation of the PCE's coefficient may not be a triv- 
ial task when the model is nonlinear and/or the dimension 
n of random input variables is relatively high. In these 
cases, the existing approaches may require a significant and 
non-trivial manipulation of the model, or the computation 
of large numbers of simulation runs, such that the advan- 
tages of the method with respect to standard MC simula- 
tions might be lost. We propose here a novel approach to 
derive the PCE's coefficients, based on convex optimiza- 
tion. Under mild assumptions, this method can be easily 
applied to any existing model, since it just requires the pre- 
liminary computation of a small number of sampled values 
of V. In this approach, a relatively large order of the poly- 
nomial chaos is initially chosen, hence a high number of 
terms in the expansion, and the PCE's coefficients are then 
computed by means of a single multi-objective optimization 
problem, exploiting £i regularization techniques (see e.g. 
Il23l . ll24l ). We provide a new, systematic way, particularly 
suited to the properties of Polynomial Chaos Expansions, 
to choose the weighting matrices in the cost function of the 
optimization problem. Moreover, we show how different 
kinds of available information on the stochastic model, in- 
cluding bounds on v and on its variance, can be easily taken 
into account as convex constraints in the optimization. As 
a result, the method is able to provide an accurate descrip- 
tion of the statistics of v{6), with few simulation runs. We 
present three case studies in a broad range of different fields, 
to demonstrate the effectiveness of the method and its ease 
of implementation. 

2 Problem Settings 

We consider a time-invariant system in state-space form: 

dt (J) 

V = v{T) = h{x{T),w{T),c,t) 

where < S M is the time variable, x{t) E M"" is the sys- 
tem state, w{t) S M"™ is an unknown input, c G W^' 
is an unknown parameter vector, finally w G M is a vari- 
able of interest, evaluated at a given time instant T. w{t) 
and c are assumed to have stochastic nature, in a sense that 



will be better detailed afterwards. Bold symbols indicate 
vectors of variables, e.g. 1 , where 

is the vector transpose operation. The aim is to derive an 
approximation of the first and second order moments and 
the pdf of V, starting from (possibly stochastic) initial con- 
ditions a;(0, c), by using the model ([U. The variable c 
accounts for uncertainty both in the model equations (e.g. 
due to uncertain physical parameters) and in the initial state 
a;(0, c), while w{t) accounts for unknown external inputs, 
like disturbances. The parameters c and the unknown in- 
put w{t), t G [0, T] in ([T]i can be typically expressed as 
functions of a n-dimensional vector 9 G K"' of indepen- 
dent and identically distributed (iid) random variables di, 
with known pdf fg, such that 0i € K d £^(0, T, P), Vi G 
{1, . . . , n}. Here, (51, P) is a probability space, is the 
set of elementary events, F is the ct— algebra of the events 
and P is the probability measure. The expectation (or first- 
order moment) of a generic random variable 6* : ^ — > M 
is denoted as i? [6*] = J^6{uj)dP{uj) = where 
Fg{k) ^ P{d < fc} is the probability distribution function 
of 9 over A. C^{il,J-, P) is the Hilbert space of all random 

1/2 

variables 9 whose La-norm, ||6l||2 ^ E , is finite, 

where | • | denotes the absolute value. K is a subspace of 
£^(f2, J^, P) that contains only centered random variables 
(i.e. ye e K, E [9] = O). Finally, the pdf of 9 is given 
by feik) — dFg / dk, and the variance (or second-order mo- 
ment) of 9 is indicated as Var(0) E [{9 ~ E[9])^] = a^, 
where ag is the standard deviation of 9. 
In many cases of practical relevance, the time-invariant pa- 
rameters c are naturally iid variables (e.g. when c stands for 
some variation of a physical parameter of the system, that 
is uncertain due to production variability). If the variables 
Ci, i G {1, Tic} have different probability distributions, it 
is possible to map a standard (i.e. with zero mean and unit 
variance) distributed Gaussian random variable, indicated 
as Af{0, 1), to a random variable with distribution func- 
tion Fc by the transformation F^^{eif{J\f{0,l))), where 
erf is the Gaussian distribution function, see e.g. 1251 . As 
regards the input w{t), t G [0, T], this can typically be 
modeled as a stochastic process or random field w{t, 9) : 
R X iiT" — >■ M"™ , which can be represented as a finite series 
of n iid random variables multiplied by deterministic func- 

n 

tions Wi{t), i G {1, n], i.e. w{t, 0) ~ ■»i'o(i) + 'u>i{t)9i 
(see for example |l26l|22l|28l[ll). 

We assume that the solution of the dynamical equations ([T]i 
in the time interval [0, T] exists and it is unique almost 
surely, i.e. with probability one. In the described context, 
the variable of interest is a random variable, v{9), and we 
name the system ([U "stochastic system". We assume that 
v{6) has finite variance: 

Assumption 1 (Finiteness of variance of v{6)) 

v{e) G c^{n,T,p). 
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Assumption[T]is typically satisfied in practical applications. 
The problem of simulating a stochastic system may be very 
complex and the main technique employed so far in engi- 
neering applications is the well-known MC approach, which 
consists in the following steps: 

Algorithm 1 (Monte Carlo simulations) 

1. extract M iid samples r € {1, M}, of 9, ac- 
cording to its distribution; 

2. for each sample, compute (or numerically simulate) 
the solution of (HJ and the corresponding value of 
i;(0M),Vre{l, M}; 

3. analyze the statistics of the collected data. 

Remark 1 For simplicity and without loss of generality, 
we consider a scalar variable v and a single value of T: 
multiple variables of interests Vj, j € {1, n„} and time 
values Ti, i £ {1, N} can be easily treated by consider- 
ing each variable and each time instant separately from the 
others, as it is done in the case studies reported in this pa- 
per. In these cases, a single simulation of the system pro- 
vides all of the corresponding samples Vj{Ti,6(^j,)), € 
{1, n4, G {1, N}. 

Although MC simulations proved to be very effective in 
many applications, the required computational times may 
be prohibitive in various cases, e.g. when a decision has to 
be taken in relatively little time on the basis of the simula- 
tions' outcome, or when repeated MC simulations have to 
be carried out to tune some input or parameter, or when the 
simulation has to be embedded in a numerical optimization 
procedure (e.g. for optimal design or control of stochastic 
systems). Polynomial Chaos Expansion techniques (see e.g. 
|3J) are able to significantly reduce the computational effort 
required by standard MC approaches, by replacing the sim- 
ulation of a (possibly very complex) dynamical system with 
the evaluation of a static function v{6) « v{6). The main 
features of PCEs are recalled in the next Section. 



3 Polynomial Chaos Expansions 

PCEs were first introduced by Wiener HI, who considered 
Gaussian random variables 9. Later on, Cameron and Mar- 
tin [2 1 showed one of the key properties of PCEs, namely 
their ability to uniformly approximate any random process 
with finite second-order moments. The polynomial chaos 
is an orthogonal basis of £^(ri, J-", P), hence any random 
variable v{9) G J^, P) has the L2-convergent expan- 

sion 111: 



;(0)=^a,cl>„,(0), 



(2) 



Table 1 : Examples of orthogonal polynomials for different 
kinds of probability measure 



Random variable 


Polynomial basis 


Gaussian 


Hermite 


Uniform 


Legendre 


Gamma 


Laguerre 


Beta 


Jacobi 



^ ■ • , E\v{e)<i>^.{e)\ 

where the coefficients a^- are given by Ofe — — (e)^]""^' 

and <i>Q,j. = $Q;fc (^'i, ■ • • , dn) is the fc-th multivariate poly- 
nomial in the series, corresponding to the /c-th vector of 
indices, or "multi-index", 0.^ = [aij^, . . . , a„^fc] , ai^k S 
N. More specifically, for a given vector of indices a^., we 

n 

have $a,(6') = Jl where $(„^^)(e'0 is the 

1=1 

univariate polynomial of degree fc, chosen according to 
the Askey scheme |4|. As an example, Hermite polynomi- 
als are used with Gaussian input random variables. Table 
[T] shows the suitable orthogonal polynomials for different 
kinds of input random variables. The choice of the univari- 
ate polynomials is made in order to satisfy the orthogonality 
property: 



(3) 



where 5ij — l\f i — j and in any other case. The coef- 
ficients of the univariate polynomials can be usually com- 
puted via a recursive equation, starting from the terms of de- 
gree and L As an example, Legendre polynomials, which 
are orthogonal w.rt. to the uniform probability distribution, 
can be obtained as: 

(4) 

n 

We denote with Ik '^i.fc the sum of the indices in 

i=l 

the multi-index a.k, and we assume that the ordering of the 
multivariate polynomials $0,^ in ©, and of the related co- 
efficients flfc, is such that Ik < ^fe+i. For practical reasons, 
the series (|2]l is truncated by considering only the multi- 
indices up to a maximal total degree I, i.e. Vctfe : Ik < I- 
An example of ordering of all the multivariate polynomials 
corresponding to 7 = 2, n = 3 is shown in Table |2] It can 
be clearly noted that the number of terms in the truncated 
series grows rapidly with n and I. Since all the possible 
multi-indices a that sum up to ^ < 7 are considered, the 
total number L of terms in the truncated expansion is: 



fc=0 



{n + l)l 



(5) 
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Table 2: Example of multivariate polynomials used in polynomial chaos, corresponding to Z = 2, n = 3 
Order Multi-index Multivariate Polynomial 








fO 01 


OLn 


(0) 


— 1 




1 




1,0,0 




.1 V / 






1 


0:2 = 


[0,1,0] 




.M 


= *(l)(^2) 




1 


"3 = 


[0,0,1] 




AO) 


-*(l)(e3) 




2 


0:4 = 


[2,0,0] 




AO) 


= *(2)(^l) 




2 


"5 = 


[0,2,0] 




AO) 


= *(2)(6'2) 




2 


ote = 


[0,0,2] 




AO) 


= *(2)(6'3) 




2 


ai = 


[1,1,0] 




AO) 


= $(l)(0l)$(l)( 


02) 


2 


as = 


[1,0,1] 




AO) 


= $(l)(0l)$(l)( 


03) 


2 




[0,1,1] 




AO) 


= $(1)(02)$(1)( 


03) 



and the series takes the form: 

L-1 

5(0) = ^afc$„,W = (6) 

where a = [ap, . . . , ol-i]"^ and 

*(0) = [$„„(0),..., (0)] (7) 

are, respectively, the vectors of the PCE's coefficients and 
of the multivariate polynomials evaluated at 6. We refer to 
the truncated expansion v{9) « v{6) ^ as the PCE of the 
random variable v{6). The PCE has been shown to con- 
verge exponentially in the L2 -sense as the maximal order 
I increases, see e.g. HJIU. By applying the orthogonal- 
ity property (|3]l, the first and second order moments of the 
random variable v{9) can be computed directly from the 
coefficients of its PCE, as follows: 

E [v{e)] = ao (8) 

a^{a) = War{v{9)) = alE (9) 

fe=i 

where ao is the coefficient of the polynomial of order I — 
(i.e. $ao = 1) in the PCE. As regards the practical 
computation of equation (|9]), the terms E [^ak (0)^] , V/c e 
{1,L ~ 1} have to be computed once for all uses, and they 
can typically be obtained quite easily. As an example, for 
Legendre polynomials and uniformly distributed input ran- 
dom variables, note that, by considering that fg — Q.5d9, 
the L2-norm squared of the multivariate polynomials is: 

n 1 

i=i 2ai,A; + 1 (10) 

Vk e {0,L- 1}. 

Similar equations can be derived for the other types of or- 
thogonal polynomials. 

Moreover, a Monte Carlo approach can be used to estimate 
the pdf of v{9) (and, hence, of v{6)) once the coefficients 
of its PCE are known, by simply evaluating the PCE v{9), 



instead of simulating the model ([T]) at step 2) of Algorithm 
[T] The computational time required to evaluate the PCE is 
often orders of magnitude smaller than the one required to 
integrate numerically the model hence the advantage of 
using polynomial chaos. 

Clearly, one of the crucial points in the use of PCEs for the 
simulation of stochastic systems is the computation of the 
expansion's coefficients, a. In the literature, this task is car- 
ried out essentially in two different ways. A first method 
(see e.g. |4|) relies on a Galerkin projection to obtain an 
augmented set of deterministic differential equations, which 
can be solved to compute the PCE coefficients. While this 
method is quite attractive from a theoretical point of view, 
it might be affected by some practical issues, since for com- 
plex nonlinear models it may be difficult and too time-con- 
suming to derive the augmented set of differential equa- 
tions, and the number of such equations may be too large 
to obtain an efficient numerical solution with standard ODE 
solvers. 

A second approach is known as Probabilistic Collocation 
Method (PCM, see e.g. ||6l |2T|), and it basically consists 
in the estimation of the coefficients from a finite number of 
data, i.e. of i' exact values of v, named "collocation points", 
corresponding to 1/ values of the input random variables, 
0(r), r e {1, I/}. Here, we consider a PCM-like approach 
for the computation of PCEs, since it appears to be more vi- 
able for the analysis of large-scale, complex stochastic dy- 
namical systems, and we propose a new method to estimate 
the coefficients. The method and its features are described 
in the next section. One of the main advantages of prob- 
abilistic collocation is that no modification to the original 
model ([T]i is required, but just a series of preliminary simu- 
lations to collect the data to be used in the coefficients' com- 
putation; one of its main disadvantages is that the number 
of collocation points can be very high, for problems with 
relatively high stochastic dimensions (i.e. high values of n) 
and strong nonlinearities. Here, we will show, through a se- 
ries of case studies, that our method yields very good results 
even with a very low number of collocation points. 
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4 Computation of Polynomial Chaos 
Expansions via convex optimization 

Given the maximal order I of the PCE and the correspond- 
ing number of terms L (|5]l, we propose the following algo- 
rithm to estimate the PCE's coefficients a: 

Algorithm 2 (PCE computation via convex optimization) 

1. sample a finite number v of independent values of the 
vector of input random variables S {1, v\, 
according to its distribution; 

2. carry out v simulations of the system ([T]), each one 
corresponding to one of the extracted samples B(r)> 

3. collect the obtained values of the variables of interest 
in the vector v = . . . , 

4. select the maximal order I for the PCE v{6) dSjl and 
compute the matrix 



where the vectors ^{6 (^rj) , r € {1, i/} , are computed 
according to (|7]l; 

5. solve the following convex optimization problem to 
compute the PCE's coefficients: 



min II Wall 1 



In (II lab , for a generic vector y e M™ the li and I2 vec- 
tor norms are defined as ||y||i = ^ \yk\ and ||y||2 — 



1/2 



X] yi I , respectively. The diagonal weighting matrix 

k=l J 

W is defined as W = diag {w{lk)) <E R^'^-^, where Ik is the 
order of the multi-index a.^, and w(Zfc), k €E {0, L — 1} is 
a sequence of scalar weights with the following properties: 

wih) > 0, Vfce {0, 
wih) > w{lj) <^ 
maxw(Zfc) = 1. 

k 



We denote with a* a global minimizer of (fTTT i and with 
v*{6) the related PCE. The convex constraints ( II Ibb are 
optional and they will be better specified later on; we now 
consider the unconstrained problem (lllab and discuss its 
features. 

Weighted £i-norin regularization of the coefficients 

The convex cost function in ( lllab is the weighted sum of 
two terms. The first one, || Wa|| 1, is a weighted ^i-norm of 
the PCE coefficients, in which the weights increase mono- 
tonically with the order of the related multivariate polyno- 
mials (see (fTzb ). £1 -norm regularization is a well-established 
technique in function approximation and regression analy- 
sis |23|,|24|, and it is a convex relaxation of the problem 
of computing an approximation which is sparsest, i.e. with 
minimal number of non-zero terms or, equivalently, with 
minimal £0 quasi-norm, see e.g. [29\. However, the use of 
the weighting matrix W is novel and pertains to the particu- 
lar properties of polynomial chaos expansions. In practice, 
minimization of the term ||Wa||i yields an estimated coef- 
ficients' vector in which the terms related to higher-order 
multivariate polynomials have smaller absolute value. The 
reason for including this term in the cost function (II lab 
is twofold: on the one hand, it accounts for the fact that, 
due to the exponential convergence property of polynomial 
chaos, the absolute values of the PCE coefficients should 
decrease as the order of the corresponding polynomials in 
the expansion increases; on the other hand, it avoids over- 
determination of the fitting problem, when the number ly of 
data is lower than the number L of coefficients. Indeed, the 
£1 -norm regularization allows one to select an initial "high" 
maximal order I, and then to rely on the convex optimiza- 
tion procedure to correctly "pick" the terms that have higher 
relevance, even if the number ly of sampled data points is 
much lower than the number of coefficients. We note that 
an approach with a similar goal, i.e. to obtain PCEs with 
low number of nonzero coefficients, has been also proposed 
in 1301 191. by means of an iterative algorithm, in which the 
maximal order I and the number v of data are gradually in- 
creased until a satisfactory accuracy is reached. Here, the 
coefficients are computed in "one shot", through the solu- 
tion of a single convex optimization problem, without any 
iterative algorithm, and the employed number of data can 
be very low. 

Weighted £2 -norm fitting of the data 

Ik > Ij, ^k,j e {0, L — 1}; xhe second term in the cost function ( lllab is the empirical 

L2-norm of the error between the PCE and the computed 



-/3||A(i-*a)||2 



subject to 
convex constraints 



(11a) 



(lib) 



k=l 



(12) 

Moreover, /? > is a scalar weight, w(-) and /3 are param- 
eters chosen by the user. Finally, the diagonal matrix A ^ 

diag contains the values A . . . , /e(^(,.))]^ 

of the pdf /e, evaluated at the considered samples 6{^r) 1 f G 
{1, 4. 



data, i.e. ||A('U - *a)||2 ~ E[\v-v\'^Y/'^ = \\v-v\\2, and 
it accounts for the convergence property, in the L2 -sense, 
of the PCE to the random variable v. The scalar weight (3 
can be used to achieve a tradeoff between the accuracy of 
the PCE, with respect to the collected data, and its complex- 
ity, in terms of weighted i'l-norm. In practical applications. 
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with a "high" value of /3 the algorithm|2]yields good results 
and it is quite robust with respect to different choices of 
weights w(^), provided that the properties (fT2b are satisfied. 

If no constraints are included, problem (fTTl l can be cast as a 
quadratic program, and a global solution can be efficiently 
computed also with thousands of coefficients f3T|. We now 
present two kinds of convex constraints ( II Ibb . which can 
be used to take into account specific additional information 
on the random variable v, at the cost of a possibly higher 
computational time. These constraints are not meant to be 
exhaustive of the possibilities that the convex optimization 
approach can open, when combined with polynomial chaos. 

Explicit maximal variance constraint 

Since the variance a"^ of the expansion is a quadratic func- 
tion of its coefficients (see eq. (|9]l), an upper bound ct^ > 
on the variance can be explicitly enforced by the convex 
constraint: 

a-^(a)<o^2. (13) 

This constraint is always feasible, since the value a — Q 
satisfies it. 

Explicit bounds on the PCE 

If some convex bounds on v are known, e.g. positiveness, 
these can be easily included in the problem as follows. As- 
sume, without loss of generality, that the bounds can be ex- 
pressed as g{v) < 0, where 5 : R — M is a convex function 
(multiple convex bounds can be reduced to this form by tak- 
ing the maximum among all of them). At first, a finite num- 
ber of further iid samples 6(^r ) of 9 has to be computed, to- 
gether with the corresponding vectors <I>(0(r)), r S {1, fj.} 
d?). Then, the following /i convex constraints can be in- 
cluded in dUbb : 



R 



5($(0(,))a) < 0, Vr e 



(14) 



Indeed, if the computed minimizer a* of ( fTTT l satisfies the 
sampled constraints (fT4l i. it is not guaranteed that the in- 
equality f{v*{9))) < is satisfied with probability one; 
however, some probabilistic results have been established in 
the context of random convex programming (see ll32l[33l ). 
and these can be used to tune the number of constraints fi. 
As a final remark, we note that, if the bounds on v are such 
that g{Q) < 0, then the convex optimization problem (fTTl l is 
feasible with probability one in the presence of the /i con- 
straints (O, since these are always satisfied by the value 
a = 0. 



5 Application examples 

In this Section, we present the results obtained by applying 
the proposed approach in three different fields. In particular, 
the first example is related to the simulation of an electric 



u(t) 



kit) 



C vc(0 



D 



Figure 1: RLC circuit with stochastic parametric uncer- 
tainty and stochastic input. Layout of the considered elec- 
tric circuit. 

circuit, with both parametric uncertainty and a stochastic 
input. The system has weak nonlinearities, it evolves in 
continuous time, it has two continuous state variables and 
13 input random variables. The second example is con- 
cerned with a model for organization innovation 134] . Such 
a model is nonlinear, it has seven positive, continuous states 
and it evolves in discrete time. The number of input random 
variables is 12. Finally, the third example is in the field of 
systems biology and presents the evaluation of the effects 
of extrinsic noise in the simulation of a chemical oscillator 
This last model is simulated through the stochastic simu- 
lation algorithm (SSA) method [351, it evolves in continu- 
ous time, it has 9 positive, discrete states and 16 input ran- 
dom variables. All together, these examples show how the 
convex optimization method can be applied in a straightfor- 
ward way to problems in a broad range of fields and with 
significant complexity, in terms of number of input random 
variables, nonUnearities, and constraints on the variables of 
interest. 

5.1 RLC circuit with stochastic parametric uncertainty 
and stochastic input 

Consider the electric circuit depicted in Fig. [T] The system 
equations are 

'idt) - -\vc{t)~^iL{t) + \u{t) 

^ (15) 

Vcit) = ■^{lL{t)-lD{t)), 

The resistance R is assumed to be a random variable R = 
iZo(l + 0.3 ^i), where Rq — 3.5 fl and 9i is a random vari- 
able with uniform distribution over [—1,1]. The inductance 
L and the capacitance C are nonlinear functions of the cur- 
rent iL{t) and voltage vc{t), respectively: 



= 0.5L_(l+exp(ai iL{t)^)) 
C{vc{t)) = 0.5 C (1 + exp(a2 vc{t)^)), 



(16) 



where ai — —0.5 10*, 02 = —0.5 10^. As an example, the 
function C{vc) is depicted in Fig. |2] Moreover, the maxi- 
mal values L and C, achieved when vc{t) = iL{t) = 0, are 
equal to L = Lo(l + 0.2 6*2), C = Co(l + 0.2 6I3), where 
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Figure 2: RLC circuit with stochastic parametric uncer- 
tainty and stochastic input. Nonlinear characteristic of ca- 
pacitance C{vc){H). 
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Figure 3: RLC circuit with stochastic parametric uncer- 
tainty and stochastic input. Comparison between the exact 
covariance function Cd {ti , ^2) of the stochastic input (sohd 
hne) and the covariance function CD{ti,t2) obtained with 
the Karhunen-Loeve expansion (dashed line). 



62, O3 are also random variables with uniform distribution 
over [—1,1]. Oi, 62, 6*3 are assumed to be independent. 
The circuit is connected to a device D, which may act both 
as load and as generator, by applying a current ioit)- In 
particular, ijj (t) is assumed to be a stochastic input of the 
form: 



■ , ^ . , 27r 

iD{t) = 0.4 sm — t 
as 



03 *D,rand(0: 

llO-^s. The 



where 03 = 1 10 ^, 04 = 5 10 ^ A, 05 
term sin ( — t] is a known sinusoidal component, while 

(t) is a random process with mean = and expo- 
nential covariance function CD{ti,t2)'- 



Cd(^i, ^2) 



(17) 



with (TD — 1 and fio — 50. In order to model the random 
process iD,rsLnd{t) as a function of a finite number of ran- 
dom input variables, we employ the Karhunen-Loeve (KL) 
expansion (see e.g. ||3]) with 10 independent random vari- 
ables 9j:i I, ... ,8jj iQ, uniformly distributed in the interval 



*D.rand(i) — io + 



10 

E 



DA 



gD,i{t)OD,i 



(18) 



In ( fTSl ), agj-, — 1/ VS is the standard deviation of the in- 
dependent random variables On.i, and Xd.i, gD.,i{t), i E 
{1, 10} are, for a given maximal time range [— T, T], the 
first ten eigenvalues and eigenfunctions of the integral equa- 
tion: 



CD{tl,t2)9DA{t2)dt2 = \DA9DA{tl)- 



(19) 



In the case of the exponential covariance function (flTl l. the 
eigenvalues Xo^i are computed as (see [SJ): 



A 



_ 2a% lip 

D.i 9 , ? ' 



where are the solutions to the following transcendental 
equations: 

— ujD,i tan(Tcj£i,i) — 0, i odd 
^D,i + tan(TcjD^i) — 0, i even. 

The corresponding eigenfunctions are: 

COs{uJD,i t) 



sm{ujD.it) 



i odd 



IT- 



sm{2T uj£)^i t) 

2 WD,: 



I , I even 



We are interested in simulating the system subject to the 
constant input u{t) = 1 10^^ V, starting from the steady 
state conditions iL{to) = OA, wl(<o) — 1 10^^ V, from 
time tg = up to time — 0.02s. Thus, we choose T — 
0.02 s to compute the KL expansion. Fig. [3] shows a com- 
parison between the exact covariance function CD{ti,t2) 
and the covariance function 6*0(^1,^2) obtained with the 

10 

KL expansion, which is computed as CdI^i, ^2) = i^D,i 

1=1 

gD,i{ti) gD.i{t2)), while Fig. |4]shows, as an example, 10 
different realizations of the signal init). We want to ana- 
lyze the statistics (first and second order moments, and pdf) 
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Figure 4: RLC circuit with stochastic parametric uncer- 
tainty and stochastic input. Examples of 10 different re- 
alizations of the stochastic input ioit) computed with the 
Karhunen-Loeve expansion. 

of the current iL{ti) and voltage vc{ti) for ti = iTg, i ^ 
{1, 10}, with Tg — 2 10^^ s. The input random variables 
6 include the 3 random variables related to parametric un- 
certainty, 9i, 62, O3, plus the 10 random variables involved 
in the KL expansion, 6*13 i £ {1, 10}. Thus, there are 
?T, = 13 input random variables in total, all uniformly dis- 
tributed in the interval [—1, 1]. According to Table [T] the 
PCE is formulated by using Legendre polynomials, whose 
coefficients can be easily computed via the recursion (|4|l. 
We applied the convex optimization procedure to estimate 
the coefficients of a different PCE for the values of vc {ti,6) 
and iL{ti, d) at all the considered time instants. In particu- 
lar, we carried out = 20 initial simulations by extracting 
the corresponding values of ^ G {1, i'} according to 
its distribution, and we used a maximal order 1 = 2 for the 
PCE. This results in L = 105 multivariate polynomials in 
the expansion. Note that the number = 20 of consid- 
ered simulations is significantly small as compared to the 
number of random input variables, n ~ 13, and to the num- 
ber of coefficients that have to be identified, L = 105. We 
made this choice to highlight one of the main advantages 
of the proposed convex optimization approach, i.e. to be 
able to obtain quite good accuracy even with an exiguous 
number of data. Standard methods, like least square fit- 
ting, do not share the same feature. In order to carry out 
a comparison between the approach proposed here and a 
standard least squares technique, we used the same data 
vc{U,9(r))i 'i'L{U,0(r)), * ^ {1, 10}, T G {1, i^} to iden- 
tify the PCE coefficients both with our convex optimiza- 
tion procedure and with the following 2-norm minimization 
problem: 

min lii; - ^alk (20) 



Problem ( |20l i is a standard least-square regression and it 
is convex, too, however in this case it is over-determined. 
Moreover, it does not take into account the available infor- 
mation on the PCE, particularly the fact that the coefficients 
related to lower-order terms are likely to be more impor- 
tant in the expansion. The PCEs obtained with the convex 
optimization approach are denoted as iL{ti,9), vc{ti,0), 
while those obtained by means of least-squares regression 
are denoted as i^'^iti, 6), v^^{ti,6), i e {1, 10}. 
In the convex optimization approach, the weights w(^), I E 
{0, 2} have been chosen as w(0) = 0.00025, w(l) = .5, 
w(2) = 1, the scalar weight /? as /? = 5 and the optimiza- 
tion problem has been solved by using the CVX package 
i36l forMatLab®. 

Fig. |5] shows the courses of the estimated mean values of 
vc{ti,6) and iL{ti,0) for all of the considered time in- 
stants U, i G {1, 10}, obtained with 100,000 MC simula- 
tions and with the PCEs tL(ti,6), vc{U,0) and i\^^{U,9), 
v^^{U,e), i e {1, 10}. We recall that, for the PCE ap- 
proximations, the first moment of the process is computed 
by simply taking, for each ij, the coefficients of the polyno- 
mial of degree in the PCEs (see ([8]l). It can be noted that 
both PCE approximations (obtained either with the convex 
optimization approach proposed here, or with least squares 
regression) give an accurate estimate of the mean of the 
variables of interest. However, the results concerning the 
variance, reported in Fig. |6] are much different: while the 
PCE identified with the convex optimization approach pro- 
posed here achieves very good results, as compared with 
the extensive MC simulations, the least squares approach 
shows a poor accuracy. In particular, the PCEs ij^^{ti, 6), 
VQ^{ti,0), i G {1, 10}, identified through the least-square 
approach, show a much lower variance with respect to the 
other two estimates. The variances of the PCEs have been 
computed with the relationship ^ (and, in particular, ( fTOb ) 
and they have been compared, in Fig. |6l to the empirical 
variance estimated by means of 100,000 MC simulations of 
the system model. 

The poor accuracy given by the PCE obtained through least 
squares regression is further highlighted by the estimates of 
the pdf, as shown as an example in Fig. |2]for the variable 
4^(^10, d) : while the empirical pdf computed with the PCE 
10 1 ^) results to be very close to the one computed with 
the standard MC approach, the one given by ij^^{tio, 0) is 
very different. Finally, as regards the computational times, 
the 100,000 MC simulations of the model ([TSll took about 
3800 s on a Intel® Core™ 2 Duo processor at 1 .3 Ghz, with 
4 GB RAM and MatLab® 2009, while the corresponding 
100,000 evaluations of the PCE required about 24 s on the 
same hardware. The average time required to solve the con- 
vex optimization problem in our procedure was 0.45 s, on 
the same computer. 
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Table 3: States, equations and parameters for the stochastic model of innovative search 



State variables 



Incoming ideas, // 
Internal stocks, IS 
New ideas, NI 
Organizational ideas, OI 
Testing ideas, TI 
Allocation of attention, AA 
External stocks, ES 



Model equations 

II{t + 1) = ci ES{t + 1) AA{t + 1) 
IS{t + 1) = C2 IS(t) + NI{t) + II{t) 
NI{t + \)^c-iIS{t) 

01{t + 1) = C6 NI{t + 1) + C4 IS{t + 1) + C5 II{t + 1) 

TI{t+l) = CTTI{t) + OI{t + l) 

AA{t + 1) = AA{t) + C8 NI{t) + eg NP{t) + cw TI{t) + cu TP{t) 
ES{t + 1) = ci2 - II{t) 



Parameter mean Parameter std. dev. 



Cl 


= 0.1375 


Cci 


= 0.0225 


C2 


= 0.2 


Cca 


= 0.02 


C3 


= 0.5 


CTC3 


= 0.06 


C4 


= 0.2 


C^C4 


= 0.02 


C5 


= 0.2 




= 0.02 


C6 


= 0.5 


CTce 


= 0.02 


C7 


= 0.275 




= 0.025 


Cs 


0.1375 


^cs 


= 0.0225 


cg 


= -0.0150 


CTcg 


= 0.002 


Cio 


= -0.0505 


O'cio 


= 0.0099 


Cll 


= 0.00055 


e'en 


= 0.00009 


Cl2 


= 1.0055 




= 0.0009 



5.2 Stochastic model of innovative search 

This second example is concerned with a dynamical model 
of how organizations pursue innovation, i.e. how they allo- 
cate attention to devise new ideas, choose part of them for 
potential investigation, and finally selects the concepts to be 
actually tested. The model has been developed by [34], it 
accounts 7 state variables and 12 stochastic parameters, and 
it has been conceived in discrete time (i.e. the state at the 
next observation time is a function of the state at the cur- 
rent observation time). An overview of the model equations 
and of the related parameters is given in Table[3] In particu- 
lar, the organization allocates a certain quantity of attention, 
AA, to the search for incoming ideas, //, from an external 
stock of ideas and information, ES. Part of these incom- 
ing ideas is stored in the internal stocks of information, IS, 
from which new ideas NI are derived. The time evolu- 
tion of the organizational ideas that are actually selected for 
possible investigation, 01, is then influenced by //, IS and 
NI. Part of 01 feeds the testing ideas, TI, i.e. the ideas 
that, among the organizational ideas, the organization actu- 
ally chooses to pursue. The state of the system is given by 



X = [//, IS, NI, OI, TI, AA, ESf and it is a vector of 
non-negative variables (i.e. x e M^, Xi > OVz G {1, 7}). 
The value of x at the observation period t + 1 is given by a 
set of six nonlinear dynamical equations and one algebraic 
equation, involving the state at the observation period t and 
12 parameters c £ R^^. The parameters are supposed to 
be independent and distributed according to Gaussian dis- 
tributions, with standard deviations and mean values 
Ci, i e {1, 12}, as shown in Table[3] Thus, in this example 
the input random variable is a vector of 12 independent 
Gaussian variables with normal distribution, such that a = 
Ci + (JciOi, i € {1, 12}. Extensive MC simulations with 
this model can be obtained with quite low computational 
times, however this example is still interesting, from the 
point of view of the approach proposed in this article, since 
the model is nonlinear and the considered parameter varia- 
tions may lead to structural changes in the stability proper- 
ties of the system (chaotic behavior), from stable modes, to 
oscillations, to divergent modes. We are interested in com- 
puting the time course, over 30 observation periods, of the 
first and second order moments and of the pdf of new ideas 
(TV/), in front of the considered variability in the model pa- 
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Figure 5: RLC circuit with stochastic parametric uncer- 
tainty and stochastic input. Mean values at t — U, i £ 
{1, 10} of (a) current iL{t) and (b) voltage vc{t) ob- 
tained either with MC simulations of the model (dashed 
line with 'o') or with the coefficients of the term of degree 
in the PCEs iL{ti, 0), vc{ti, 0) (solid line with '*') and 
i^^iU, 6), Vc^iU, 6) (dotted line with '[>')■ 

rameters. As a matter of fact, for the observation periods 
t = 1, 2 and 3 the value of NI remains fixed at its initial 
condition, so only the periods t = 4, . . . , 30 are analyzed. 
Following our convex optimization procedure, we sample 
f — 300 values of 6 E [—1,1]^^ and compute the cor- 
responding values of NI{t, 0), t E {4, 30}, starting from 
the same initial condition a;(0) = [0, 0, 0, 0, 0.2, 50]^. 
We use Hermite polynomials (according to Table[T]) and we 
consider a maximal order 1=3 for the PCEs NI{t, 9), 
so that each expansion has L — 455 terms. We note that a 
Galerkin projection method would lead, in this case, to 3 1 85 
discrete-time dynamical equations (i.e. seven model equa- 
tions, times 455 coefficients in the PCE), while a standard 
least-squares regression would need at least 455 sampled 
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Figure 6: RLC circuit with stochastic parametric uncer- 
tainty and stochastic input. Variances at < ^ ti, i E {1, 10} 
of (a) current iL{t) and (b) voltage vc{t) obtained either 
with MC simulations of the model (dashed line with 'o') or 
with the coefficients of the PCEs iL{ti, 6), vdU, 0) (solid 
line with '*') and i'l^{ti,6), v^^{ti,6) (dotted line with 
'O'). 

data, to avoid over-determination. We select the weights 
w(0, / E {0, 3} as w(0) = 0.0001, w(/) = ^V^ E {1, 3} 
and the scalar weight (3 — 10^. Moreover, since variable 
NI is defined to be positive, we include 500 additional con- 
straints NI(t,6(^r)) > 0, Vr e {1, 500}, as described in 
the main paper. Finally, we include a quadratic constraint 
on the expansions' variances, by setting ct^ = 2 a^, where 
(T is the empiric standard deviation computed by using the 
300 simulated data. We solve the convex optimization prob- 
lem again with the CVX package. The courses of the first 
and second order moments of NI{t, 6), computed by using 
the PCEs' coefficients via (l8]l-(|9]l, are shown in Fig. |8ja)- 
(b), where they are compared with the empirical moments 
obtained by means of 100,000 MC simulations. Fig. |9] 
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Current (A) ^ jq-5 Current (A) ^ Current (A) ^ ^^-5 



Figure 7: RLC circuit with stochastic parametric uncer- 
tainty and stochastic input. Comparison between the em- 
pirical pdfs of variable iL{tw,S) estimated with 100,000 
MC simulations with the system model (left) and with 
100,000 MC evaluations of the PCEs iL{tio,0) (middle) 
and li^iho,0) (right). 

shows the time evolution of the pdf of NI, estimated either 
by computing 100,000 MC simulations or with the corre- 
sponding 100,000 evaluations of the PCEs. Finally, Table|4] 
shows, as an example, a comparison between the quartiles 
of the variables of interest, computed either with standard 
MC simulations, or with the corresponding PCE. The good 
matching between all of these statistics shows that, also in 
this case, the expansions' coefficients, computed with the 
proposed convex optimization method, are able to describe 
the stochastic process with good accuracy.This second ex- 
ample is concerned with a dynamical model of how organi- 
zations pursue innovation, i.e. how they allocate attention 
to devise new ideas, choose part of them for potential in- 
vestigation, and finally selects the concepts to be actually 
tested. The model has been developed by |34|, it accounts 
7 state variables and 12 stochastic parameters, and it has 
been conceived in discrete time (i.e. the state at the next 
observation time is a function of the state at the current ob- 
servation time). An overview of the model equations and 
of the related parameters is given in Table |3] In particu- 
lar, the organization allocates a certain quantity of attention, 
AA, to the search for incoming ideas, //, from an external 
stock of ideas and information, ES. Part of these incom- 
ing ideas is stored in the internal stocks of information, IS, 
from which new ideas NI are derived. The time evolu- 
tion of the organizational ideas that are actually selected for 
possible investigation, 01, is then influenced by //, IS and 
NI. Part of 01 feeds the testing ideas, TI, i.e. the ideas 
that, among the organizational ideas, the organization actu- 
ally chooses to pursue. The state of the system is given by 
X = [II, IS, NI, 01, TI, AA, ESf and it is a vector of 
non-negative variables (i.e. x G M7 , Xi > OVi G {1, 7}). 
The value of x at the observation period t + lis given by a 



Table 4: Stochastic model of innovative search: quartiles 
of the probability distribution of new ideas during 27 ob- 
servation periods, estimated either with standard MC sim- 
ulations, or with the PCEs computed with the convex opti- 
mization approach 

Obs. MC simulations Polynomial chaos 
period 25% 50% 75% 25% 50% 75% 



set of six nonlinear dynamical equations and one algebraic 
equation, involving the state at the observation period t and 
12 parameters c € M}'^. The parameters are supposed to 
be independent and distributed according to Gaussian dis- 
tributions, with standard deviations and mean values 
c,;, i e {1, 12}, as shown in Table[3] Thus, in this example 
the input random variable is a vector of 12 independent 
Gaussian variables with normal distribution, such that Ci — 
Ci + ac-di, i E {1, 12}. Extensive MC simulations with 
this model can be obtained with quite low computational 
times, however this example is still interesting, from the 
point of view of the approach proposed in this article, since 
the model is nonlinear and the considered parameter varia- 
tions may lead to structural changes in the stability proper- 
ties of the system (chaotic behavior), from stable modes, to 
oscillations, to divergent modes. We are interested in com- 
puting the time course, over 30 observation periods, of the 
first and second order moments and of the pdf of new ideas 
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(NI), in front of the considered variability in the model pa- 
rameters. As a matter of fact, for the observation periods 
t = 1, 2 and 3 the value of NI remains fixed at its initial 
condition, so only the periods t — 4, . . . , 30 are analyzed. 
Following our convex optimization procedure, we sample 
ly ~ 300 values of E [—1,1]^^ and compute the cor- 
responding values of NI{t, 6), t E {4, 30}, starting from 
the same initial condition a;(0) = [0, 0, 0, 0, 0.2, 50]^. 
We use Hermite polynomials (according to Table[T]) and we 
consider a maximal order 7=3 for the PCEs NI{t, 9), 
so that each expansion has L — 455 terms. We note that a 

(a) 
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Figure 8: Stochastic model of innovative search, (a) Mean 
values and (b) variances of the New Ideas NI at the obser- 
vation periods t = A, . . . , 30, estimated either with 100,000 
standard MC simulations (dashed line with 'o') or with the 
polynomial chaos expansions computed with the convex 
optimization approach (solid lines with '*'). 

Galerkin projection method would lead, in this case, to 3 1 85 
discrete-time dynamical equations (i.e. seven model equa- 
tions, times 455 coefficients in the PCE), while a standard 
least-squares regression would need at least 455 sampled 
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New Ideas 

Figure 9: Stochastic model of innovative search. Level 
curves of the pdf of new ideas NI as a function of the obser- 
vation period, estimated by means of (a) 100,000 standard 
MC simulations and (b) polynomial chaos expansions com- 
puted with the convex optimization approach. 

data, to avoid over-determination. We select the weights 
w(0, / e {0, 3} as w(0) = 0.0001, w(/) = ^V^ e {1, 3} 
and the scalar weight /3 — 10^. Moreover, since variable 
NI is defined to be positive, we include 500 additional con- 
straints NI(t,9(^r)) > 0, Vr e {1, 500}, as described in 
the main paper. Finally, we include a quadratic constraint 
on the expansions' variances, by setting ct^ = 2 a^, where 
(J is the empiric standard deviation computed by using the 
300 simulated data. We solve the convex optimization prob- 
lem again with the CVX package. The courses of the first 
and second order moments of NI{t, 0), computed by using 
the PCEs' coefficients via (IH]!-®, are shown in Fig. ISa)- 
(b), where they are compared with the empirical moments 
obtained by means of 100,000 MC simulations. Fig. |9] 
shows the time evolution of the pdf of NI, estimated either 
by computing 100,000 MC simulations or with the coiTe- 
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spending 100,000 evaluations of the PCEs. Finally, Table|4] 
shows, as an example, a comparison between the quartiles 
of the variables of interest, computed either with standard 
MC simulations, or with the corresponding PCE. The good 
matching between all of these statistics shows that, also in 
this case, the expansions' coefficients, computed with the 
proposed convex optimization method, are able to describe 
the stochastic process with good accuracy. 

5.3 Chemical oscillator 

The last application example is in the field of systems biol- 
ogy. We consider the chemical oscillator analyzed in [371 
and we simulate this system by means of Gillespie's stochas- 
tic simulation algorithm ll35l (SSA) and the common reac- 
tion path (CRP) method proposed in ||38l . In this chemi- 
cal process, the promoters Pr and Pa control a repressor 
protein R and an activator protein A, respectively. The lat- 
ter is able to combine with either or Pa, giving rise 
to the complexes Pha and Paa, to enhance the transcrip- 
tion of mRNAA or mRNAj^, respectively, with the con- 
sequent synthesis of new A or R proteins. The repres- 
sor protein R is able to combine with A, by forming the 
intermediate complex An, and to induce its degradation. 
The state of this model is given by the quantities of the in- 
volved molecules, i.e. X = [A, R, Pa, Pr, Paa, Pra, 
mRNAA, mRNAR, Ar]'^ e IR>o, which are discrete by 
definition, while the model evolves in continuous time. In 
particular, according to the SSA-CRP simulation method, 
each one of the 16 reactions that may occur in this process 
has its own internal clock, and its own stream of random 
firing times, whose total number, at a given time instant, 
is a Poisson random variable. We indicate these streams of 
random firing times as ^ — {^i = {T'fci}fe°=0' * ^ l^}}, 
where Tk- is the time interval between two subsequent fir- 
ing times for the reaction, and fc^ e N is a counter giv- 
ing the total number of reactions of type i that already took 
place. The internal clocks of the reactions evolve at dif- 
ferent speeds, given by the propensities a^, i = 1, . . . , 16 
times the common time variable, t. The latter are gener- 
ally nonlinear functions of the state and of 16 model pa- 
rameters Ci, i — 1, . . . , 16. The model's chemical reac- 
tions, the related parameters and the propensities' equa- 
tions are described in Table |5] As an example, consider 
the 4* reaction and assume it took place already 10 times, 

i.e. fc4 = 10. When the 4* internal clock hits its own next 

11 

firing time, given by ^ r^^, the 4* reaction takes place 

fc4=0 

again, the counter ki is augmented by one, and the sys- 
tem state is updated according to the corresponding chem- 
ical equation (i.e. the number of mRNAR molecules is 
augmented by one), as well as the values of the propensi- 
ties. Then, the simulation continues with the new propen- 
sity values (i.e. the new clocks' "speeds"). Since each re- 



action, when it takes place, influences the propensities of 
the other reactions, the simulation must be carried out in 
a sequential fashion, by iteratively computing the reaction 
that fires next and updating the state, propensities and coun- 
ters (for more details, the interested reader is referred to 
e.g. [38]). Typically, SSA simulations are carried out for 
a given initial state and fixed parameters' values, by taking 
multiple random extractions of the firing times' streams ^ 
(internal noise), and then some statistic of interest is ana- 
lyzed. Yet, the model parameters are not fixed and known, 
rather they can be assumed to be random variables them- 
selves, the so-called extrinsic noise, and it is of interest to 
study the sensitivity of the SSA outcome to such parameter 
variations. PCEs have been already applied in the context 
of systems biology IfTSl . by using a projection method and 
then a quadrature approach to identify the PCE coefficients. 
As it is also recalled in |18|, Gauss-Hermite quadrature is 
efficient up to 3-5 stochastic dimensions. In this example, 
we analyze the sensitivity to random perturbations in all of 
the 16 parameters, thus we have 16 stochastic dimensions. 
Each one of the 16 model parameters Ci, i — 1, . . . , 16, is 
assumed to be uniformly distributed in the interval ±10% 
centered at the corresponding mean value as indicated 
in Table |5] Therefore, the input random variables are given 
by a vector 6 = {Oi S [-1, 1] : q = (1 + 0.16*^)^, i = 
1, . . . , 16} of 16 independent variables, each one with uni- 
form distribution in the interval [—1,1]. The variables of 
interest are the expected value A{t, 9) — E[A{t, 6, ^)] and 
the variance a\{t,e) = E[{A{t,e,£,) - ^{t,e,C)Y] of the 
number of molecules of protein A, evaluated every 5 s up to 
t = 50 s of simulation time, starting from the initial condi- 
tion x{Q) = [0, 177, 1, 1, 0, 0, 4, 0, 279]J. More specifi- 
cally, we consider the empirical values of A{t, 6), a\{t, 9), 
computed by averaging over 1000 SSA simulations. By us- 
ing the CPR method ||38l , each SSA realization is associ- 
ated to its own, fixed seed that generates the random streams 
^. In this way, each simulation is evaluated with different 
values of 9 (i.e. extrinsic noise) but always with the same 
random stream ^ (i.e. internal noise): therefore, extrinsic 
noise and internal noise are effectively decoupled, since in 
practice, for a fixed value of 6, the process is totally deter- 
ministic and given by the SSA simulations, each one with 
its own stream of random firing times. Thus, the only source 
of randomness lies in the model parameters c G R^^. In- 
deed, the application of Galerkin projection methods ap- 
pears to be not trivial in this case, due to the particularity 
of the described SSA method and to the discrete nature of 
the state variables. Moreover, we consider PCEs of or- 
der 3, which, in a 16-dimensional space, involve 969 terms. 
With the convex optimization approach, we run 100 sets of 
1,000 SSA simulations, corresponding to 100 samples of 
the random parameter vector This number of data is very 
low with respect to the dimensionality of the random pa- 
rameter 9 E [—1, 1]^^, yet the resulting 3'^''-order PCE is 



13 



Table 5: Reactions, propensities and nominal parameter values for the chemical oscillator The stochastic model parame- 
ters are uniformly distributed in the interval ±10% around the nominal value 



Reaction 


Propensity 


Nominal 
parameters 


Pa-^Pa + tuRNAa 


ai = 


ci Pa 


ci = 50 


Paa ^ Paa + tuRNAa 


a2 = 


ci5 ci Paa 


C2 = 0.01 


Pr^Pr + m.RNAR 


as = 


C2Pr 


C3 = 50 


Pra ^ Pra + mRNAn 


di = 


Cl6 C2 -TRA 


C4 = 5 


tuRNAa tuRNAa + A 


as = 


cs tuRNAa 


C5 =20 


tuRNAr tuRNAr + R 


ae = 


C4 tuRNAr 


C6 = 1 


A + R^ A- R 


ay = 


c^AR 


C7 = 50 


Pa + A^ Paa 


as = 


cqPaA 


Cg = 1 


Paa ^Pa + A 


ag = 


ctPaa 


eg = 100 


Pr + A^Pra 


aio = 


--c^PrA 


Clo = 1 


Pra^Pr + A 


an = 


-- cg Pra 


cii = 0.2 


A^0 


ai2 = 


-- cioA 


Cl2 = 10 


R^0 


ai3 = 


= cii R 


ci3 = 0.5 


tuRNAa ^ 


ai4 = 


= Ci2 rni?iVA^ 


Ci4 = 1 


tuRNAr ^ 


Ol5 = 


= Ci3 tuRNAr 


Cl5 = 10 


Ar > 


ai6 = 


- Cl4 


ci6 = 5, 000 



Table 6: Chemical oscillator: quartiles of the probability distribution of the expected number A{t, 6) of molecules of 
protein A, computed every 5 s up to 50 s, estimated either with standard MC simulations, or with the PCEs computed with 
the convex optimization approach 





MC simulations 


Polynomial chaos 


Time (s) 


25% 


50% 


75% 


25% 


50% 


75% 


5 


58 


113 


217 


51 


123 


242 


10 


547 


686 


833 


556 


693 


841 


15 


258 


325 


392 


251 


319 


388 


20 


25 


38 


56 


29 


44 


61 


25 


3 


5 


8 


4 


7 


9 


30 


16 


30 


54 


17 


33 


58 


35 


174 


263 


381 


181 


272 


386 


40 


403 


489 


591 


404 


493 


595 


45 


200 


280 


371 


196 


278 


373 


50 


37 


63 


104 


42 


67 


110 
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Table 7: Chemical oscillator: quartiles of the probability distribution of the variance (j\{t, 6) of the number of molecules 
of protein A, computed every 5 s up to 50 s, estimated either with standard MC simulations, or with the PCEs computed 
with the convex optimization approach 

MC simulations (x 10'') Polynomial chaos (x IQ-^) 



Time Ts) 


25% 


50% 


75% 


25% 


50% 


75% 


5 


54.4 


115.2 


229.3 


53.7 


129.2 


259.5 


10 


298.5 


364.7 


436.4 


294.4 


356.5 


422.6 


15 


192.7 


235.6 


285.8 


187.2 


238.5 


295.9 


20 


24.3 


37.4 


53.5 


28.4 


42.0 


58.7 


25 


3.0 


5.2 


7.8 


3.7 


5.7 


8.1 


30 


17.1 


32.9 


62.3 


19.3 


38.1 


68.9 


35 


159.3 


238.8 


346.0 


161.1 


242.5 


344.2 


40 


264.3 


335.9 


424.5 


264.1 


335.3 


422.1 


45 


164.7 


226.4 


305.0 


164.5 


231.1 


311.2 


50 


35.6 


61.1 


103.2 


41.4 


66.4 


109.2 



highly accurate with respect to the results obtained with a 
standard MC approach, considering 10,000 sets of 1,000 
SSA simulations. In particular, the comparison between the 
two approaches is shown in Fig. [T0lfT3l and in Tables |6l|7] 
In the convex optimization problem, we chose the weights 
w(0, I e {0, 3} as w(0) = 0.0001, w(/) = ^yi e 
{1, 3}, and the scalar weight (3 — 10'^. We solved the con- 
vex optimization problem by using the Yalmip ||39l tool- 
box for MatLab®. Moreover, similarly to the second ap- 
plication example, we included 5,000 additional constraints 

'I{tJ(^r)),a-l{t,e,^r)) > 0, Vr e {1, 5, 000}, in order to 
take into account the positiveness of the variables of inter- 
est. As regards the computational times, the time required 
to compute the 100 data points used to identify the PCEs' 
coefficient was 67 min, the solution of the 20 convex op- 
timization problems (2 variables of interest evaluated at 10 
different time instants) took 2 hours (averagely 6 minutes 
for each PCE), finally the evaluation of 10,000 MC values 
of the resulting PCE took 8 s on an Intel® Core^'^ 2 Duo 
processor at 1.3 GHz, with 4 GB RAM and MatLab® 2009. 
Thus, the PCE-convex optimization method took about 3 
hours in total. The time required to compute the 10000 
standard MC simulations was about 6670 minutes, i.e. 4.6 
days. The model equations for the SSA simulations have 
been programmed in Simulink®, and the computation have 
been carried out on a Speedgoat® real-time machine, by 
using Embedded Matlab® and xPC-target® tools to auto- 
matically generate the simulation code from the Simulink 
model. Indeed, the results of this example confirm that the 
proposed method, based on convex optimization, is able to 
compute the PCE's coefficients with good accuracy also in 
the presence of a relatively large number of random dimen- 
sions and model nonlinearities, with a very limited num- 
ber of preliminary data. We note that, once the PCE's co- 
efficients have been computed, 100,000 evaluations of the 
expansion would take, on the same Intel® Core^'^ 2 Duo 



computer, about 80 s, while the corresponding simulations 
with the dedicated real-time hardware would take about 46 
days. 

6 Conclusions 

We proposed a new method to compute polynomial chaos 
expansions, by means of a suitably defined convex opti- 
mization problem. The method can easily handle thousands 
of terms in the PCE, corresponding for example to stochas- 
tic dimensions of 15-20 with orders of 3-4. Bounds on the 
first and second order moments and on the values of the 
resulting PCE can also be explicitly included. We applied 
the approach to three examples in a broad range of different 
fields: in all cases, the derived PCEs, computed via a very 
low number of preliminary simulations, accurately captured 
the process' statistics, despite the presence of nonlinearities 
and high stochastic dimensions. This aspect indicates that a 
quite small number of sampled simulations already contains 
sufficient information on the process, to derive an accurate 
PCE approximation. This method can be straightforwardly 
used in a large variety of applications, since it does not re- 
quire any modification to the existing model, but just a small 
number of simulation runs. 
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Figure 10: Chemical oscillator, (a) Mean values and (b) 
variances of the expected number A{t, 6) of molecules of 
protein A estimated every 5 s up to 50 s of simulation, either 
with 10,000 standard MC simulations (dashed line with 'o') 
or with the polynomial chaos expansions computed with the 
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Figure 11: Chemical oscillator (a) Mean values and (b) 
variances of the variance a\ {t, G) of molecules of protein 
A estimated every 5 s up to 50 s of simulation, either with 
10,000 standard MC simulations (dashed line with 'o') or 
with the polynomial chaos expansions computed with the 
convex optimization approach (solid lines with '*'). 
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Figure 12: Chemical oscillator. Level curves of the pdf 
of the average value A{t,6) of the number of protein A 
molecules as a function of time, estimated by means of 
(a) 10,000 standard MC simulations and (b) polynomial 
chaos expansions computed with the convex optimization 
approach. 
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Figure 13: Chemical oscillator Level curves of the pdf of 
the variance (T\{t, 6) of the number of protein A molecules 
as a function of time, estimated by means of (a) 10,000 stan- 
dard MC simulations and (b) polynomial chaos expansions 
computed with the convex optimization approach. 
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